home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / src / demos / audio / ameshC / power.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-02  |  6.6 KB  |  277 lines
  1. /*
  2.  * Copyright 1991, 1992, 1993, 1994, Silicon Graphics, Inc.
  3.  * All Rights Reserved.
  4.  *
  5.  * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
  6.  * the contents of this file may not be disclosed to third parties, copied or
  7.  * duplicated in any form, in whole or in part, without the prior written
  8.  * permission of Silicon Graphics, Inc.
  9.  *
  10.  * RESTRICTED RIGHTS LEGEND:
  11.  * Use, duplication or disclosure by the Government is subject to restrictions
  12.  * as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
  13.  * and Computer Software clause at DFARS 252.227-7013, and/or in similar or
  14.  * successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
  15.  * rights reserved under the Copyright Laws of the United States.
  16.  */
  17. /*
  18.  * Copyright (C) 1991, Silicon Graphics, Inc.
  19.  * All Rights Reserved.
  20.  */
  21. #include <math.h>
  22. #include "fft.h"
  23. #include <malloc.h>
  24. #include <stdio.h>
  25. /*
  26. ** Instructions:
  27. **
  28. **     First call buildtable to create the sin/cos tables.
  29. **        buildtable(ln)  where ln is the log base 2 of the number of
  30. **                        samples.
  31. **     Next call bitreverse to bitreverse the input data.
  32. **        bitreverse(x,ln) where ln is the log base 2 of the number of
  33. **                         samples and x is a complex array which holds
  34. **                         the samples.
  35. **     Finally, call fft to perform a forward fft, or call ifft to
  36. **                         perform an inverse fft.
  37. **        fft(x,ln)
  38. **        ifft(x,ln)       where ln is the log base 2 of the number of
  39. **                         samples and x is a complex array which holds
  40. **                         the samples.
  41. **     If you are performing several fft/iffts of the same size, you only
  42. **                         need to build the tables once
  43. */
  44. #include        <math.h>
  45.  
  46.  
  47. /*******
  48. //  Object
  49. // To read in a sampled sound file 
  50. // compute an estimate of the power Spectrum of a Delta t
  51. // write out the pse for each time interval
  52. // pse should be normalized to 0 to 1.0
  53. // Question should we plot log(power) or power???
  54. /*******/
  55.  
  56. static complex *c=0;
  57.  
  58. #define PI2 6.283185307179586476925287
  59.  
  60. int pwr_two[] = {
  61. 1,   2,   4,   8,
  62. 16,   32,   64,   128,
  63. 256,   512,   1024,   2048,
  64. 4096,   8192,   16384,   32768,
  65. 65536,   131072,   262144,   524288,
  66. 1048576,   2097152,   4194304,   8388608,
  67. 16777216,   33554432,   67108864,   134217728,
  68. 268435456,   536870912,   1073741824
  69. };
  70.  
  71. void NRICfft( float *, unsigned int, int);
  72.  
  73.  
  74.  
  75. #define FILTERSIZE  16384
  76. #define FILTERSIZE2 32768
  77. #define LOGFILTERSIZE2 15
  78.  
  79. #define PI2 6.283185307179586476925287
  80. int nfreqs, nintervals;
  81.  
  82. #define MAXFREQS 256
  83. #define MAXINTERVALS 2000
  84. #define MAXSAMPS 8192
  85.  
  86. complex  tmpc[8192];
  87.  
  88. typedef signed short SAMPLE;
  89.  
  90. int fftinitialized =0;
  91.  
  92.  
  93. float spectrum( SAMPLE * samples, int nsamples, float * ps, int nfreqs, float maxpwr )
  94.     {
  95.     int i,lognsamples,nfft;
  96.     SAMPLE * pb;
  97.     float pwr,nfftsqd;
  98.     complex * pc;
  99.     nfft = 1;
  100.     lognsamples = 0;
  101.         while (nfft < nsamples) { nfft<<=1; lognsamples++;}
  102.     if ( fftinitialized != lognsamples)
  103.         {
  104.         buildtable(lognsamples);
  105.         fftinitialized = lognsamples;
  106.         }
  107.  
  108.  
  109.     for(pb = samples, pc = tmpc, i = 0; i < nsamples ; i+=1)
  110.         { pc->r  = (float) *pb++;
  111.           pc->i = 0.0; pc++;
  112.         } 
  113.     if (nfft > nsamples )   /* zero fill /***/
  114.        for(i = 0; i < nfft-nsamples ; i++)
  115.         { pc->r = 0.0; pc->i = 0.0 ; pc++; }
  116.     pc = tmpc;
  117.  
  118.     bitreverse(tmpc,lognsamples);
  119.     fft(tmpc,lognsamples);
  120.      nfftsqd = (float)(nfft)*(float)(nfft);
  121.  
  122.     for(i=1; i<nfreqs; i++)
  123.         {
  124.         
  125.         pwr = tmpc[i].r * tmpc[i].r;
  126.         pwr += tmpc[i].i * tmpc[i].i;
  127.         pwr = fsqrt(pwr/nfftsqd) ;
  128.         ps[i-1] = pwr;
  129.         if (pwr > maxpwr) maxpwr = pwr;
  130.         }
  131.     return(maxpwr);
  132.     }
  133.  
  134.  
  135.  
  136.  
  137.  
  138.  
  139. void buildtable(int ln)
  140. {
  141.         int n,halfn,i;
  142.         float theta;
  143.         n=pwr_two[ln];
  144.         halfn=pwr_two[ln-1];
  145.         if(c)free(c);
  146.         c=(complex *)malloc(halfn*sizeof(complex));
  147.         if(c==0){printf("malloc failed\n");exit(0);}
  148.         for(i=0;i<halfn;i++) {
  149.                 theta= PI2*i/n;
  150.                 c[i].r=cos(theta);
  151.                 c[i].i= -sin(theta);
  152.         }
  153. }
  154.  
  155. void bitreverse(complex *x,int ln)
  156. {
  157.     int i,j,k,n;
  158.     complex hold;
  159.  
  160.     n=pwr_two[ln];
  161.     k = 0;
  162.     j=0;
  163.     for (;;) {
  164.         if (k > j) {
  165.             hold.r = x[j].r;
  166.             hold.i = x[j].i;
  167.  
  168.             x[j].r = x[k].r;
  169.             x[j].i = x[k].i;
  170.  
  171.             x[k].r = hold.r;
  172.             x[k].i = hold.i;
  173.         }
  174.         j++;
  175.         if(j==n)break;
  176.  
  177.         /* k = k+1, (k is n bits long and bit reversed). */
  178.  
  179.         i=ln-1;
  180.         while(i>=0) {
  181.             if (k & pwr_two[i])
  182.                 k = k - pwr_two[i];
  183.             else
  184.                 break;
  185.             i--;
  186.         }
  187.  
  188.         k = k + pwr_two[i];
  189.     }
  190. }
  191.  
  192. void fft(complex *x, int ln)
  193. {
  194.         unsigned int stuckbit,stuckbitcomplement,halfn,t,b;
  195.         unsigned int shft,msk,wt;
  196.         unsigned int c0,c1;
  197.         unsigned int n;
  198.         float qr, qi, xtr, xti, xbr, xbi, cwtr, cwti;
  199.  
  200.         n=pwr_two[ln];
  201.         msk=n-1;
  202.         halfn=pwr_two[ln-1];
  203.         shft=ln-1;
  204.  
  205.         stuckbit=1;
  206.         c0=ln;
  207.         while(c0--){
  208.             t=0;
  209.             stuckbitcomplement = ~stuckbit;
  210.             c1=halfn;
  211.             while(c1--) {
  212.                 wt=(t<<shft)&msk;
  213.                 b=t+stuckbit;
  214.  
  215.                 xbr = x[b].r;
  216.                 xbi = x[b].i;
  217.                 cwtr = c[wt].r;
  218.                 cwti = c[wt].i;
  219.  
  220.                 qr = xbr * cwtr - xbi * cwti;
  221.                 qi = xbr * cwti + xbi * cwtr;
  222.  
  223.                 xtr = x[t].r;
  224.                 xti = x[t].i;
  225.  
  226.                 x[b].r = xtr-qr;
  227.                 x[b].i = xti-qi;
  228.  
  229.                 x[t].r = xtr + qr;
  230.                 x[t].i = xti + qi;
  231.  
  232.                 t=(b+1) & stuckbitcomplement;
  233.             }
  234.             stuckbit <<= 1;
  235.             shft -= 1;
  236.         }
  237. }
  238.  
  239. void ifft(complex *x,int ln)
  240. {
  241.         unsigned int stuckbit,stuckbitcomplement,halfn,t,b;
  242.         unsigned int shft,msk,wt;
  243.         unsigned int c0,c1;
  244.         unsigned int n;
  245.         float qr,qi;
  246.  
  247.         n=pwr_two[ln];
  248.         msk=n-1;
  249.         halfn=pwr_two[ln-1];
  250.         shft=ln-1;
  251.  
  252.         stuckbit=1;
  253.         c0=ln;
  254.         while(c0--){
  255.             t=0;
  256.             stuckbitcomplement = ~stuckbit;
  257.             c1=halfn;
  258.             while(c1--){
  259.                 b=t+stuckbit;
  260.                 wt=(t<<shft)&msk;
  261.  
  262.                 qr =     x[b].r * c[wt].r  + x[b].i * c[wt].i;
  263.                 qi =   -(x[b].r * c[wt].i) + x[b].i * c[wt].r;
  264.  
  265.                 x[b].r = x[t].r-qr;
  266.                 x[b].i = x[t].i-qi;
  267.  
  268.                 x[t].r += qr;
  269.                 x[t].i += qi;
  270.  
  271.                 t=(b+1) & stuckbitcomplement;
  272.             }
  273.             stuckbit <<= 1;
  274.             shft -= 1;
  275.         }
  276. }
  277.